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The Rayleigh capillary instability of a cylindrical interface between two immiscible 
fluids is one of the most fundamental in fluid dynamics. As Plateau observed from energetic 
considerations and Rayleigh clarified through hydrodynamics, such an interface is linearly 
unstable to fission due to surface tension. In traditional descriptions of this instability 
it occurs everywhere along the cylinder at once, triggered by infinitesimal perturbations. 
Here we explore in detail a recently conjectured alternate scenario for this instability: front 
propagation. Using boundary integral techniques for Stokes flow, we provide numerical 
evidence that the viscous Rayleigh instability can indeed spread behind a front moving at 
constant velocity, in some cases leading to a periodic sequence of pinching events. These 
basic results are in quantitative agreement with the marginal stability criterion, yet there 
are important qualitative differences associated with the discontinuous nature of droplet 
fission. A number of experiments immediately suggest themselves in light of these results. 
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Recently, Bar-Ziv and Moses discovered the "pearling instability" of lipid bilayer mem- 
brane tubes, in which the application of laser tweezers to the membrane induces a periodic 
modulation in the radius [|l|. Theoretical explanations for this behavior invoke a ten- 
sion created in the membrane by the laser tweezers. Were this an interface between two 
immiscible fluids, such tension would induce capillary breakup of the cylinder into droplets 
via the Rayleigh instability In the pearling phenomenon however, breakup is prevented 
by the membrane bending elasticity. A more striking difference is that the pearling in- 
stability propagates — the modulated state is observed to invade the uniform cylindrical 
region at a constant velocity. This propagation is totally unlike traditional descriptions of 
the Rayleigh instability , in which perturbations grow uniformly along the length of the 
tube. 

In a previous publication , we conjectured that the viscous Rayleigh instability can 
propagate as well. Propagating fronts in which a new stable state invades an unstable 
region are a common feature of overdamped systems , arising in such diverse situations 
as reaction-diffusion systems, dendritic growth, and, more recently, type-I superconductors 
. These previously well-studied examples all display a continuous evolution, in contrast 
to the discontinuous evolution associated with drop fission. Below we show by numerical 
methods that despite this fundamental difference, the Rayleigh instability of a cylindrical 
interface can indeed propagate. Furthermore, the coarse features of this behavior can be 
described by the marginal stability criterion (MSG), an analytical method which has found 
wide applicability in the simpler situations mentioned above [§. 

The discontinuous evolution associated with rupture of the thread can have significant 
implications for the possible existence of a propagating front. As a droplet pinches off 
from the main body, the two tips of the broken neck recede from the pinching point; if 
the retracting end overtakes the front, propagation will be spoiled. Experiments on the 
breakup and relaxation of elongated drops suspended in an outer fluid reveal that this 
competition depends on the relative viscosity of the two fluids. When the drop viscosity 
r]~ is much smaller than the outer viscosity 77+ , the drop breaks before its ends have time 
to retract; in the other extreme, the ends retract significantly before breaking off (for 
long enough drops). We shall see below that when the viscosity ratio X = r]~ /r]~^ is such 
that retraction is slow on the time scale for breakup, there is a propagating front moving 
with constant velocity (Fig. 1). No analytical solution is known for the complex shape 
evolution shown in Fig. 1. Indeed, this is the case in most examples of front propagation, 
where, however, the MSG nearly always provides correct predictions for front speed. In 
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our studies of drop breakup below we find that the front speed is rather accurately given 
by the linear MSG, a remarkable result in light of the strongly nonlinear, singular shape 
evolution behind the front. In addition to the front velocity, we compute the time between 
primary pinching events in the breakup of the cylinder using the MSG quantities and find 
good agreement with the numerical calculations. 
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Fig. 1: Sequence of drop shapes for viscosity ratio A = 0.05 at = 
6.67n?7+i?/7; n = 1,2,3, ...,15 from top to bottom. To illustrate the 
complete evolution, we have drawn the daughter droplets (but not the 
satellite droplets). However, the evolution of each connected component 
was computed independently of the others. 
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Drop evolution. We look for the propagation of the Rayleigh instability in a numer- 
ical simulation of a long, axisymmetric, approximately cylindrical drop, initially stationary 
and with hemispherical caps on the ends. The axis of the cylinder defines the x-direction, 
so that the shape is given by the surface of revolution generated by the radius r{x) We 
work in the overdamped limit in which the inertial terms of the Navier-Stokes equations 
may be disregarded. Therefore, the outer and inner fluids with velocities and pressures 
are described by the Stokes equations and the constraint of incompressibility 

r^±V^u± = Vp^, V-u± = 0. (1) 

Propagation in the inertial regime is complicated by the presence of dispersive capillary 
waves, and we leave this case and the case of net flow (like a jet) for future work. The 
boundary conditions at the interface 5" are continuity of fluid velocity u+l^ = u~|g, con- 
tinuity of tangential stress, and the jump in normal stress: 

^.(4-^^^)ls = -27^^n„ (2) 

where n is the outward surface normal, H is the mean curvature, 7 is the (constant) 
interfacial tension, and the stress tensors are = Tj^iy + V juf) — P^^ij- For an 
axisymmetric shape with radius r(x), and with = dr/dx, the mean curvature is 
H = (r^^/2)(l + r^)~"^/^ — (2r)~^(l -|- r^)"^*^^. The only important material parameters 
are therefore the surface tension 7 and the viscosities r]~^ and 77". 

Note that (|^) gives rise to the absolute (rather than convective) instability of a sta- 
tionary cylindrical interface. When the fluid velocity is zero, the pressure jump across the 
interface is Ap = 2'yH. Consider an axisymmetric perturbation in which the radius is 
slightly pinched to form a neck. If the disturbance is of sufficiently long wavelength, the 
magnitude of the curvature is increased, thus raising the pressure in the neck. Fluid is 
therefore forced out of the pinched region, leading to growth of the perturbation. 

The challenging numerical task of solving the three-dimensional Stokes equations is 
greatly simplified by the boundary integral technique, in which (|ip and (|^) are recast as an 
integral equation for quantities on the interface []T^. This approach removes one spatial 
dimension from the problem, leading to 
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+ ^i^P.V.^t.,(x)n,(x)K,^.,(x,y)d^(x), (3) 



where y is a point on the interface S, u(y) is the velocity of the interface at y, V.V. denotes 



the principal value, and the kernels J^j and K-j^ are Green's functions [|TT |. 

Axisymmetry reduces our problem to an integral equation with one space and one 
time dimension. This equation is solved numerically by standard adaptive-grid techniques 
||12|| to yield the interfacial velocities, and so the drop shape, as a function of time. This 
procedure can describe the continuous motion of the interface only; to describe the breakup 
of the interface we simply demand that pinching occurs whenever the radius r(x) < 0.005R, 
where R is the initial cylinder radius. Although somewhat arbitrary, we expect this choice 
for the cutoff leads to little error in the gross evolution since further decreases in the neck 
radius occur rapidly owing to the large velocity and curvature gradients near the pinch 
point. Furthermore, the size of the region in which these quantities grow large is small 
In any case, these same factors make it difficult to do numerical calculations when the neck 
approaches rupture. Once a droplet has pinched off, we neglect it in further calculations, 
since experience has shown the the evolution of the droplet that pinches off has little effect 
on the evolution of the main drop fl^ . 



Our findings for the shape evolution are in qualitative accord with previous investi- 
gations 0. At higher values of A, i.e. around A = 10, the evolution is dominated by 
retraction. As an example. Fig. 2 shows that the A = 10 drop simply retracts and shows 
no sign of developing a neck in the time it takes the A = 0.1 drop to break a few times. 
Fig. 3 shows a magnified view of the evolution of these droplets. Note that at a given 
time the radius r{x) of the A = 10 drop grows monotonically in x to its maximum value, 
but the radius for A = 0.1 is modulated; these undulations are the seeds of the capillary 
instability. 

When the outer viscosity is not too small, we can estimate the retraction speed by 
assuming the bulge on the drop end is spherical and balancing tension with drag 0. 
This leads to a velocity that scales like ■y/r]'^, since the drag is controlled by the outer 
viscosity [^. This estimate only accounts for the dissipation due to the flows in the two 
fluids induced by dragging a spherical interface through the outer fluid; it disregards the 
contribution due to shape change, which become comparable when A > 1. Indeed, when 
there is no outer fluid all the dissipation takes place inside the retracting end. The estimate 
is borne out by our simulations: if we take the retraction velocity to be the speed of the 
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Fig. 2: Numerical results. Propagation of the Rayleigh instability for 
A = 0.1 (left). Retraction-dominated dynamics for A = 10. (right). 

= 6.Q7nr]~^R/j, where n — 1,2,3, ...,15 labels the shapes from the 
top. 

maximum of r{x) nearest the end of the drop (Fig. 4), we find retraction velocities that 
range from 0.27/77+ at A = 0.005 to 0.17/77+ for A = 10. 

The front region connects the uniform cylinder, which has not yet undergone the 
Rayleigh instability, to the growing bulge, which eventually pinches off (Fig. 4). We extract 
the front from the sequence of drop shapes by choosing for each shape a window with one 
edge such that r{x) is within numerical accuracy of the unperturbed radius R, and the 

5 



other edge with |r(x)—i2|/-R reaching a small value a* (typically 0.1-0.15). The drop radius 
in this window is fit [|16| with an exponential envelope times a sinusoidal function of the 
form r{x) = 0^^ + 03 exp{—q"x) cos{q'x — a^). We find that the fit parameters a^^, Og, q', q" 
do not change with small changes in the choice of fitting window or small changes in the 
initial shape of the ends of the drop. The front speed is the velocity of the point Xj^{t) for 
which the envelope function 03 exp(— .(t)) = a*. 




Fig. 3: Enlarged view of the droplet ends of Fig. 2 for A = 0.1 and 
A = 10. = (37.5 + 3n)r]^R/'y, where n = 1, 2, 6 labels the shapes. 

Turning now to the analysis of the fronts, we find a smoothly moving front for a 
range of viscosity ratios, from A = 0.005 to A = 1.0. For example, the top graph in Fig. 5 
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shows the front position as a function of time for A = 0.05. At low values of A, around 
A = 0.005 |T^, the ends do not retract much before the drop pinches (Fig. 2). This is the 
"end-pinching" behavior described in ||^ : the daughter droplets pinch off one at a time in 
a periodic fashion at the ends of the main drop [0. Despite the discrete nature of these 
events, there is a front moving smoothly at constant velocity. We therefore interpret this 
end-pinching at low viscosity ratio as a front of the Rayleigh instability. 

As A increases, retraction becomes more apparent in the time it takes the droplets 
to break off. When A reaches a value at which the retraction speed becomes comparable 
to the speed at which the pinches are spreading (around A = 1), we find that the drop 
evolution can not be described as a smoothly moving front. The details of this process are 
complicated. For A ~ 0.2, the periodic behavior of the breakup ceases. The droplets that 
pinch off are not of uniform size and do not pinch off at a constant rate. Nevertheless, we 
are still able to fit the advancing profile to a front moving smoothly with constant velocity. 
When A > 1.0, this smoothly propagating front behavior is lost. Around this value of A we 
also observe necks that start to form at one point but then heal, leading to breakup at a 
different point along the drop. This behavior is reminiscent of that seen in the experiments 
of Tjahjadi et a/., who observed that a drop stretched to a cylindrical shape just beyond 
the critical aspect ratio for breakup does not break fl^. Instead, two necks form and then 
disappear, and the elongated drop retracts to a sphere without breaking. We can offer 
no explanation of this behavior, and merely note that when retraction is faster than the 
pinching process, there is no reason to expect the simple picture of the Rayleigh instability 
outlined in the last section to apply. 

Front velocity and the MSG. For completeness, we review first the basic facts of 
the marginal stability criterion. In the leading edge of the front, the amplitude w(x, t) = 
r{x,t) — R is small and can be represented as a linear combination of Fourier modes each 
having the form u{x, t) ~ exp{u!{q*)t + iq'*x), in which u{q), possibly complex, is the linear 
growth rate, and the real and imaginary parts of the wavevector q = q' + iq" describe the 
periodicity and sharpness of the mode. The velocity of the exponential envelope of any 
mode is = Reiv/lmq. The MSG selects a unique mode q* (and hence velocity v*) by the 
condition that the envelopes of modes nearby in q neither outrun nor fall behind that of 
q* . This provides two additional relationships between v and g, leading to the conditions 

V* = Im— - =0, V* = Re^, (4) 

Img* oq oq 
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Fig. 4: A snapshot of the end of the main drop for A = 0.05, in which 
the aspect ratio has been distorted to show the front more clearly. The 
center of the drop is at x = 0. The darkened line is the front region; 
the inset is a magnified view of this region (solid line) and a fit to the 
front (dotted line). The retraction speed is taken to be the velocity of 
the maximum marked with a plus sign. 



where uo* = uj{q*) A rigorous proof of the validity of this approach is known 

only for the Fisher-Kolmogorov equation = 



+ u 



for which it can be 



shown that sufficiently localized initial conditions will evolve into a front moving at a 
unique speed v = v* [^. A very large number of partial differential equations have been 



shown to be correctly described by this principle or its nonlinear variants 0, even when 
no analytic front solution is known. There is, however, no general criterion to determine 
when the MSG is applicable, and thus the question of its validity in this new setting of 
discontinuous evolution is of great interest. 

The MSG quantities depend only on the growth rate uj, which is known from To- 
motika's generalization ||2^ of Rayleigh's result to the two-fluid problem. For a cylin- 
der of radius R, uj is a, function of wavevector q and viscosity ratio A in the form 
uj{q,X) = {■y / Rr]~^)A{q, X){1 — q^) . The dynamical factor A{q,X), too lengthy to quote 
here, accounts for viscous dissipation, and the factor 7(1 — g^) is associated with the 
Young-Laplace force due to a constant- volume distortion of wave number q. Inserting this 
growth rate into the MSG equations (^) yields the front velocity previously derived in 
as a function of viscosity ratio. Fig. 6 compares this prediction with our present numerical 
calculations. 

Note that since the growth rate for the Rayleigh instability does not contain the 
physical process of retraction, we expect the MSG velocity to depart from the true front 
speed when the front and retraction speed are comparable. This may be the source of 
disagreement between the MSG and the numerical simulations for A > 0.5. As A is increased 
beyond this value, retraction becomes more important until flnally the simple picture of 
a propagating Rayleigh instability becomes invalid. At the lower end of the viscosity 
contrasts we studied, the boundary integral method starts to lose accuracy; since it is 
difficult to quantify these errors we cannot say if this is the source of disagreement between 
the MSG and the numerical calculations ||24|| . 

While the front position advances smoothly in accord with the MSG picture, the 
behavior of the front shape encoded in the wavenumber q' and the inverse front width q" 
is more complicated. We find that these quantities are not constant but rather oscillate 
with approximately the period of pinching; however, they are always near the MSG values. 
We show one example in the bottom graph of Fig. 5. We conclude that the MSG captures 
the gross features of the front shape but not the flne structure. 

Given that we flnd the breakup process to be periodic in the viscosity range A = 0.005- 
0.1, it is natural to ask if the MSG can predict the time tpi^ch between primary pinching 
events, once the transients due to the initial shape settle down. A plausible but naive 
starting point is to use the shortest characteristic growth time, tpinch = '^(?max)~^' where 
9max wavenumber of the fastest growing mode. However, it is well known that the 

fastest growing mode is not the mode that is selected when the MSG is applicable. In 
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Fig. 5: Top: results of boundary integral calculations and fitting routine 
for front position vs. time. Bottom: comparison of prediction of MSG 
and results of boundary integral method for the initial front wavenumber 
q' and front width q" . Both plots are for A = 0.05. The simulation begins 
at t = 0; by ^7/(77+/?) = 20 transients from the initial shape have died 
out. 

the standard MSG treatment, there are in fact two selected wavenumbers: q', the selected 
wavenumber in the leading edgeof the front, and = lm{uj{q*) +iq* v*) / v* [^, the selected 
wavenumber in the saturated pattern. Thus another candidate for the pinching time is 
^pinch — 27^/(^0^*)- However, as we pointed out in [||, we do not expect to be relevant 
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Fig. 6: Front velocity v vs. A. The circles are the results of the numerical 
calculations, the line is the prediction of the MSG from [§. 

to the Rayleigh problem since the saturated pattern has not continuously evolved from the 
leading edge of the front. Therefore we expect the appropriate wavelength to be 27: /q'. In 
Fig. 7 we show that the numerical simulations clearly favor the relation tpi^ch = 2n/{q'v*). 

In conclusion, we have shown that the breakup of elongated drops can proceed via 
front propagation for the viscosity ratios in the range A = 0.005-1.0. Remarkably, despite 
the complex nature of the discontinuous dynamics, the marginal stability criterion gives 
a good description of the gross quantitative features of the shape evolution, such as the 
front velocity and shape, and the pinch-off time. There are many experimental systems 
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Fig. 7: Various predictions for the primary pinch time. Sohd hne is 
2Tx/{q'v), dotted hne is 2Tx/{qQv), and the dashed hne is l/uj{q^^^. The 
circles are from the simulation. The breakup is aperiodic for A > 0.2, 
leading to a distribution of droplet sizes. 

which could be used to test our results. The most natural candidate is the four-roll mill, 
in which a planar hyperbolic flow stretches a spherical drop into a long thread 0. When 
this flow ceases, the relaxation or breakup of the thread can be studied. Other possibilities 
are capillary bridges stabilized by electric or magnetic flelds [25| or very viscous jets. 
Experimental conflrmation of our picture of propagating topological transitions would be 
an important contribution to the understanding of interface motion in fluid dynamics. 
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